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Abstract 

We  consider  the  use  of  preconditioning  methods  to  accelerate  the  convergence 
to  a  steady  state  foi  both  the  incompressible  and  compressible  fluid  dynamic  eipia- 
tions.  Most  of  the  analysis  relies  on  the  inviscid  ecpiations  though  some  ai)plications 
for  viscous  flow  are  considered.  The  preconditioning  can  consist  of  either  a  matrix 
or  a  differential  operator  acting  on  the  time  derivatives.  Hence,  in  the  steady  state 
the  original  steady  solution  is  obtained.  For  finite  difference  methods  the  precondi¬ 
tioning  can  change  and  improve  the  steady  state  solutions.  Several  preconditioners 
previously  discussed  are  reviewed  and  some  new  a})])roarhe.s  are  presented. 


‘  I'liis  research  wa.s  supported  l)y  the  National  Aeronautics  and  .Space  Administration  under  NA.SA 
Contract  No.  NASl-lOdSO  while  the  author  was  in  residence  at  the  Institute  for  Computer  Applications 
in  Science  and  Engineering  (ICASE),  NASA  bangley  Research  Center,  Hampton,  V'A  23681-0001. 
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1  Introduction 


Over  tlie  past  years  luiineroiis  researchers  have  tried  to  solve  the  steady  state  iiieonipress- 
ihle  ecpiatioiis  for  both  iiiviseid  and  viscous  flows.  This  also  lead  to  attempts  to  solve  the 
compr('ssil)h'  ('(luations  over  a  large  range  of  mach  numbers.  A  standard  way  of  solving 
the  steady  stat('  equations  is  to  march  the  time  dependent  equations  until  a  steady  state 
is  reached.  Since  the  transient  is  not  of  any  interest  one  can  use  acceleration  technicpies 
which  destroy  the  time  accuracy  but  enable  one  to  reach  the  steady  state  faster.  For  the 
incompresible  erpiations  the  continuity  equation  does  not  contain  any  time  derivatives. 
To  overcome  this  difficulty  (fliorin  [I  t]  added  an  artificial  time  derivative  of  the  pressure 
to  the  continuity  e(p.iation  together  with  a  multiplicative  variable,  /:?  ,  With  this  artifi¬ 
cial  term  the  resultant  scheme  is  a  symmetric  hyperbolic  system  for  the  in  viscid  terms, 
"fhus,  th('  system  is  wi'll  poserl  and  and  numerical  method  for  hyperbolic  systems  can 
Ix'  used  to  advance  this  system  in  time..  The  free  parameter  ji  is  then  chosen  to  reach 
the  steady  stati'  cpiickly.  I.ater  Turkel  [54]  ('xtended  this  conci'pt  by  adding  the  pressure 
time  derivative  to  the  momentum  equations  and  introducing  a  second  free  parameter  a. 
"fhis  systc'm  can  then  be  analyzed  for  optimal  o,  /T  The  resulting  system  after  precon¬ 
ditioning  is  no  longer  symmetric  but  can  be  symmetrized  by  a  change  of  variables.  This 
will  be  shown  in  more  detail  later. 

It  is  well  known  that  it  is  difficult  to  solve  the  compressible  equations  for  low  Mach 
numbers.  For  an  explicit  schenie  this  is  easily  seen  by  looking  at  the  time  steps.  For 
stability  the  time  step  must  be  chosen  inversely  proportional  to  the  largest  eigenvalue  of 
th('  system  which  is  approximately  the  speed  of  sound,  c,  for  slow  flows.  However,  other 
waves  are  convected  at  the  fluid  sp<*etl,  u  ,  which  is  much  slower.  Hence,  these  waves 
<lon't  change  very  much  over  a  time  step.  Thus,  thousands  of  time  steps  are  required 
to  rc'ach  a  steady  state.  Should  one  try  a  multigrid  acceleration  one  finds  that  the  same 
disparity  in  wave  speeds  slows  down  the  multigrid  acceh'ration.  With  an  implicit  method 
an  ADI  factorization  is  usually  used  so  that  one  can  easily  invert  the  im])licit  factors. 
The  use  of  ADI  introduces  factorization  errors  which  again  slows  down  the  convergence 
rat(’  when  there  are  wave  speeds  of  very  different  magnitudes  [49]  . 

F(U’  small  Mach  nundjers  it  can  be  shown  ([28],  [41]  )  that  the  incompressibleequations 
approximate  the  compressible  equations.  Hence,  one  needs  to  justify  the  use  of  the 
compressible  equations  for  low  Mach  flows.  We  present  .several  reasons  one  would  still 
use  the  compressible  equations  even  though  the  Mach  number  of  the  flow  is  small. 

•  Tlu're  are  many  sophisticated  compressible  codes  available  that  could  be  used  for 
such  jjroblems  especially  in  complicated  geometries 

•  For  low  speed  aerodynamic  problems  at  high  angle  of  attack  most  of  the  of  the 
flow  consists  of  a  low  Mach  number  flow.  However,  there  are  localized  regions 
containing  shocks. 

•  In  many  problems  thermal  effects  are  important  and  the  energy  equation  is  coupled 
to  the  other  ('quations. 

Therefore',  one  wants  to  change  the  transient  nature  of  the  system  to  remove  this 
disparity  of  the  wave  speeds.  Ffased  on  an  analogy  with  conjugate'  graelie'iit  methe)els 
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smli  inetliods  were  called  [54]  preconditioned  methods  since  the  object  is  to  reduce  the 
condition  number  of  the  matrix.  Another  approach,  in  one  dimension,  is  to  diagonalize 
the  matrix  of  the  inviscid  term.  One  can  then  use  a  different  time  stej)  for  (-ach  equation, 
or  wav<'.  Upon  returning  to  the  original  variables  one  finds  that  this  is  etpiivah'iit  to 
multiplying  the  time  derivatives  by  a  matrix.  Hence,  this  same  api)roach  is  named 
characteristic  time  stepping  in  [55].  In  multidimensions  one  can  no  longer  comph'tely 
decouple  th('  waves  by  diagonalizing  both  the  entropy  and  the  shear  waves  and  so  the 
characteristic  time  sti'pping  is  oidy  an  approximation. 

Thus,  for  both  the  incompressible  and  compre.ssible  equations  we  will  consider  sys- 
t(>ms  of  the  form 

( 1 )  +  fj-  +  =0, 

This  system  is  written  in  conservation  though  for  some  applications  this  is  not  necessary. 
Onr  analysis  will  be  based  on  the  linearized  ecpiations  so  that  the  conservation  form  does 
not  appear  in  the  analysis  though  it  does  appear  in  the  numerical  system.  This  syteui 
is  now  replaced  by 

(-)  P  '  +  /)■  +  <fy  =  0. 

or  in  linearized  form 

(■?)  P~'  ((',  -f  Ati'j.  +  Bwy  =  0, 

In  order  for  this  system  to  be  eqiuvalent  to  the  original  system  in  tlu'  steady  state 
we  demand  that  P  hav<’  an  inverse.  This  only  need  be  true  in  the  flow  regime  under 
consifleration.  We  shall  see  later  that  frequently  P  is  singular  at  stagnation  points  and 
also  along  the  sonic  line.  Thus,  we  will  only  consider  strictly  sid)sonic  flow  without  a 
stagnation  point  or  else  strictly  supersonic  flow.  For  tran.sonic  flow  it  is  necessary  to 
smooth  out  the  singularity  in  a  neighborhood  eff  the  sonic  line.  We  also  assume  that  the 
.lacobian  matrices  A  =  ^  and  ^  are  simultaneously  symmetrizable.  In  terms  of 

the  “symmetrizing'  variables  we  also  demand  that  P  be  positive  definite.  We  shall  show 
later  in  (h’tail  that  it  does  not  matter  which  set  of  dependent  variables  are  used  to  develop 
tin'  preconditioner.  One  can  transform  between  any  two  sets,  of  variables.  The  choice 
of  variables  is  dictated  only  by  convenience  in  constructing  the  preconditioner.  Popular 
choices  are  two  out  of  density,  pre.ssure.  enthalpy,  entropy  or  temperature  in  addition  to 
the  velocity  components.  Thus,  when  we  are  finished  we  will  analyze  a  system  which  is 
similar  to  (5)  where  the  matrices  A  and  14  are  symmetri<’  and  P  is  both  symmetric  and 
positive  definite.  Such  .systems  are  known  as  symmetric  hyperbolic  systems.  One  can 
tlu'ii  multit)ly  this  system  by  w  and  integrate  by  parts  to  get  estimates  for  the  integral  of 
irf.  i.e.  energy  estimates.  These  <\stimate.s  can  then  be  used  to  show  that  tlie  system  is 
well  posed.  We  stress  that  if  P  is  not  positive  then  we  change'  the  physics  of  the  probh-m. 
For  ('xample,  if  P  =  —I  then  we’  have’  rever.se’d  the  time’  elire’ction  anel  must  the’ie’fore’ 
e  hange’  all  the-  benuielay  e'onditie)ns.  Hence,  to  be  sure  that  the  syste’ui  is  we’ll  pe).se’el  with 
the-  original  type-  eef  beeundary  e'e)nelit ieens  we  shall  eudy  e'eensieler  the’  .symuu’trie'  hype-rbeelic 
sysle’Ui.  For  meere’  ge’iie’ral  syste’ins  eeiie  must  use  a  more’  cejuijdicaleei  analysis  te)  show 
we'lb|)ose’elne’ss  feer  the’  init ial-bounelary  value’  predih’iu  ([50],  [fid]). 


With  these  assumptions  we  see  that  the  steady  state  solutions  of  tlie  two  systems 
are  the  same.  Assuming  the  steady  state  has  a  unique  solution  it  does  not  matter  which 
system  we  march  to  a  steady  state.  We  shall  later  see  that  for  the  finite  difference 
a])proximations  the  steady  state  solutions  are  not  the  necessarily  same  and  usually  the 
preconditioned  system  leads  to  a  better  behaved  steady  state. 

We  ran  also  look  at  (.'3)  from  a  different  viewpoint.  We  assume  that  the  matrices  A 
and  B  are  symmetric  and  P  is  positive  definite.  It  is  well  known  that  for  the  Euler  equa¬ 
tions  that  the  matrices  A  and  B  cannot  be  simultaneously  diagonalized  by  a  similarity 
transformation.  However,  the  matrix  P  has  changed  the  equation.  Since  P  is  positive 
definite  there  exists  a  matrix  Q  so  that  P  =  QQ*.  We  then  introduce  a  new  variable 
u’  =  Qx\  For  constant  coefficients  A,  B  (3)  is  replaced  by 

( 1)  -f  Q*AQvx  +  Q*  BQvy  —  0, 

Thus,  the  diagonalization  question  changes  and  we  wish  to  know  if  A  and  B  can 
be  simultaneously  diagonalized  by  a  congruence  transformation  {Q*AQ)  .  A  sufficient 
condition  for  this  to  be  true  is  that  there  exist  numbers  so  that  ui\A  A  ^028  is 

positive  definite.  It  is  shown  in  [53]  that  this  true  for  su])er.sonic  flow.  Hence,  we  have 
shown  that  for  supersonic  flow  one  can  introduce  a  preconditioning  matrix  so  that  the 
equations  (constant  coefficients)  are  diagonalized.  However,  this  is  not  true  for  subsonic 
flow.  We  shall  later  show  that  using  differential  operators  one  can  diagonalize  the  system 
even  for  subsonic  flow. 

2  Incompressible  equations 

We  first  consider  the  incompressible  inviscid  equations  in  i)rimitive  variables. 

Ux  +  Vy  =  0 

(5)  Ut  +  mix  +  vuy  A  px  =  0 

Vt  A  nvx  +  e?’y  +  py  =  0 

We  consider  generalizations  of  (’horin’s  pseudo-compressibility  method  [Mj.  Using  the 
preconditioning  suggested  in  [54]  we  have 


'^Pi  +  fix  +  —  0 

A 

nil 

(fi)  —xit  A  mix  A  I’Uy  A  Px  =  0 

A 

nv 

-—vi  A  uvx  +  vvy  A  Py  =  0 
A 

or  in  conservation  form 


I 

—^Pi  A  Ux  A  Vy  = 


3 


0 


(7) 


0 


(a  +  l)u  ,  ,  2  .  N  ,  /  X 

- — - ut  +  {u  +p)x  +  {uv)y  = 

(a  +  1  ,  ,2 

- +  {tiv)^  +  (v^  +  p)y  =  0 

Hence,  (7)  reduces  to  the  original  pseudo-compressibility  method  when  a  =  0.  The 
conservative  form  reduces  to  the  basic  method  when  cv  =  —  1  .  We  can  also  write  (7)  in 
matrix  form  using 

/  1  0  0  \  /  fP  0  0  \ 

(8)  Pt~^  =  I  oiuj/3^  1  0  1  Px  =  I  —au  1  0 

\  avjfP  0  1/  \  —av  0  1  y 

i.e. 

/0  10\/p\  /00l\/p 

-fluoj  uj-|-0u0  u 

^^0  0  u  j  \1  0  11  )  \  V 

Multiplying  by  P  we  rewrite  this  as 

(10)  wt  +  PAwx  +  P  Bwy  —  0 

We  also  define 

(11)  D  —  U>\  A  Ll?2  B  —  1  ^  ^1 X  ^2  —  1 

where  u!\,uj-2  are  the  Fourier  transform  variables  in  the  x  and  y  directions  respectively. 
The  speeds  of  the  waves  are  now  governed  by  the  roots  of  det{XI  —  PAu)\  —  PBijj2)  —  0 
or  equivalently  de/-(AP“'  —  Aui  —  Buj2)  —  0.  Let 

(12)  <7  =  uu>i  +  vuz 
Then  the  eigenvalues  of  D  are 

(13)  do  =  <7 

d±  =  \/2  ( 1  -  a)<7  ±  ^/(l  -  a^q'^  -f  4/3^ 

Note  that  in  the  special  case  a  =  1  we  have 

(14)  d±  =  ±/^ 

and  so  the  ‘acoustic’  speed  is  isotropic. 

We  see  that  the  spatial  derivatives  involve  symmetric  matrices,  i.e.  D  is  a  symmetric 
matrix.  Thus,  while  the  original  system  was  symmetric  hyperbolic  the  preconditioned 
system  is  no  longer  symmetric.  In  ([■'>4])  it  is  shown  that  as  long  as 

(15)  >  o{n^  +  v^) 


(  l//3^  0  0  \  /  p 

(9)  j  au/d^  10  j  u 
\  av/d^  0  1  /  \  V 
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then  the  system  is  symmetrizable.  Hence,  for  any  nonnegative  a  the  system  is  always 
symmetrizable.  Recall  that  o  =  0  for  the  original  pseudo-compressibility  equations  in 
primitive  variables  (7)  while  a  =  — 1  for  the  original  pseudo-compressibility  method  in 
conservative  variables  (8)  For  o  =  1  we  need 


(16) 


IP  >  [u^  +  v^) 


On  the  other  hand  the  eigenvalues  are  most  equalized  if  iP  =  (u^  +  v^).  Hence,  we  wish 
to  choose  (P  slightly  larger  than  -f  v^.  However,  numerous  calculations  verify  that  in 
general  a  constant  (3  is  the  best  for  the  convergence  rate.  The  reasons  for  this  are  not 
clear. 

However,  we  wish  to  stress  that  (3  has  the  dimensions  of  a  speed.  Therefore,  [3  can  not 
be  a  universal  constant.  There  are  papers  that  claim  that  (3  =  I  or  f3  =  2.5  are  optimal. 
Such  claims  can  not  be  true  in  general.  It  is  simple  to  see  that  if  one  nondimensionalizes 
the  equation  then  (3  gets  divided  by  a  reference  velocity.  Hence,  the  optimal  ‘constant’  (3 
depends  on  the  dimensionalization  of  the  problem  and  in  particular  depends  on  the  inflow 
conditions.  In  most  calculations  the  inflow  mass  is  fixed  at  one  or  else  p+{u^  +  v^)/2  =  1. 
Such  conditions  will  give  an  optimal  f3  close  to  one.  However,  if  one  chose  the  incoming 
mass  as  ten  then  the  optimal  0  would  be  closer  to  ten. 

Van  Leer,  Lee  and  Roe  considered  the  compressible  equations.  They  wanted  a  sym¬ 
metric  preconditioner  so  that  there  would  be  no  question  of  well  posedness.  We  now 
translate  their  results  to  the  incompressible  equations  (1).  They  assume  that  the  flow  is 
aligned  with  the  x  direction  and  so  v=0  and  |u|^  is  the  total  speed  of  the  fluid.  Their 
preconditioner  in  this  coordinate  system  is 


(17) 


P  = 


-fzU 
\  0 


0  \ 

1  +  ^  0 
0  r  / 


Choosing  r  =  1  preserves  the  speed  of  the  shear  wave  while  choosing  0  =  1  gives  an 
isotropic  ‘acoustic’  wave  (20)  the  magnitude  of  this  acoustic  wave  is  determined.  In 
order  to  compare  this  formula  with  the  previous  formulas  we  wish  to  reformulate  this 
preconditioner  for  the  case  where  the  flow  is  not  aligned  in  the  x  direction.  We  denote 
the  matrices  in  the  streamwise  and  perpendicular  directions  as  A\\  and  zfi  respectively. 
We  next  define  the  rotation  matrices  as 


( 1 

0 

( 

1 

0 

0  \ 

0 

cosO 

sin9 

IJ-'  = 

0 

cos  9 

—sin9 

VO 

—sinO  cosO  j 

\ 

0 

sui9 

cos9  j 

To  get  the  streamwise  direction  we  shall  choose 


cosO 


u 


sinO  — 


V 

\/  iP  +  tP 


One  can  then  verify  that  given  the  original  matrices  A,  B. 


(18)  /l|l  =  U{Acose  +  Bsine)iJ-^ 

Ax.  =  I'i—Asind-j-  Bcos6)U~^ 

Given  numbers  Lb\,i^2  for  we  define 


—  u-xsinO 
u-i  —  U\sin6  +  Cj-xcosO 


note 


Also  define 


Then  it  is  easy  to  verify  that 


“f" 


2 

1 


—  tJj  +  ui^  ■ 


p  =  f/-’p(/. 


P(A||t:;,  +  A^Qjx)  =  V  (P(Au;,  +  B^x)] 


Therefore,  the  appropriate  preconditioner  is  P  given  by 


(19) 


Pv 


^  IL^  +  —U  —V  '' 

“  *  +  u2+„2  u2+,;2 

UV  1  I 

\  n^+u^  '  ) 


Note  that  P,A,B  are  symmetric  matrices.  This  does  not  imply  that  PA  or  PB  are 
symmetric.  However,  this  is  still  a  symmetric  hyperbolic  system  and  so  the  standard 
energy  estimates  prove  the  well  posedness  of  the  system.  We  also  see  that  the  eigen¬ 
values  do  not  change  if  we  use  the  streamwise  direction  or  the  full  2D  form.  Thus,  the 
eigenvalue.s  of  the  preconditioned  system  are 


(20)  du  —  \/xB  -P  =  ucui  -t-  ruJ2  =  7 

d±  —  -p  id 

(i±  are  the  same  as  in  (13)  if  we  choo.se  o  =  1  and  3  —  \/iB  +  xd. 

As  noted  before,  with  the  preconditioner  of  Van  Leer  et.  al.  one  cannot  have  tl',e 
usual  shear  speed  together  with  an  i.sotropic  ‘acoustic’  wave  speed  with  an  arbitrary  mag¬ 
nitude.  With  therefore,  consider  a  modification  of  their  preconditioner.  In  streamwise 
coordinates  it  is  given  by 


21 


with 


P  = 


V  0 


-A  0  \ 
a  0 
0  1  / 


.  .B  ±  dyJiP  -  X?  .  2A 
a  = - ,  a  -  — , 

XI  ll 


u  0 


6 


pre 


Choosing  (5^  —  gives  the  original  preconditioner  of  Van  Leer  et.  al.  for  inconi- 
ssible  flow.  In  general  nonaligned  coordinates  this  becomes 


(22) 


Pv 


M 


( 

—au  1  + 
—av 


—ou 

(2or  — 1 
(2nt  — 1 


a  = 


IP  ±  I^^IP  -  {u^  +  v^) 


—av 

(2a— l)mt 

I  I 


+  7^  0 


tP  -h 

Now,  we  have  the  condition  fP  >  +  tP  (cf.  16).  The  speeds  are  now  given  by 

do  —  q 


d±  =  ±fi 

This  can  now  be  compared  with  (14)  for  Pj. 

NumeroTis  computer  runs  have  shown  that  P-j  works  best  with  (3  constant  and  not 
depending  on  the  speed.  To  date  there  have  been  no  computer  calculations  for  the 
incompressible  equations  with  Py- 

These  examples  show  that  the  preconditioning  is  not  unique.  If  fact,  it  is  straightfor¬ 
ward  to  see  that  the  transpose  of  Px  is  also  a  preconditioner  with  the  same  eigenvalues 
for  the  preconditioned  system.  In  general,  these  various  systems  will  have  similar  eigen¬ 
values  but  different  eigenvectors  for  the  preconditioned  system.  Numerous  calculations 
show  that  the  system  given  by  Pj  is  more  robust  and  converges  faster  than  that  with 
the  transpose  preconditioner.  This  shows  that  it  is  not  sufficient  to  consider  just  the 
eigenvalues  but  somehow  the  eigenvectors  are  also  of  importance. 


3  Compressible  equations 

The  time  dependent  Euler  equations  can  be  written  as 


(23) 
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The  first  general  attempt  to  replace  this  by  other  systems  of  equations  with  the  same 
steady  state  was  by  Viviand  ([59], [27]).  He  considered  both  incompressible  and  com¬ 
pressible  isoenthalpic  flow.  VVe  will  consider  preconditionings  that  are  a  generalization 
of  (9) 
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Note  that  if  we  use  ^  instead  of  dp  tlie  matrices  become  symmetric.  We  next  present 
the  eigenvalues  of  P  D  (defined  in  (11).  Let 

(21)  q  =  i/u;]  +  vijj-i 

then 

(25)  do  =  q 

f4  =  1/2  ( 1  -  ft  +  li'^lc^)q  ±  y/((  1  -  ft  +  iP/c^Yq^  +4(1-  r/Vc^/i'^ 

If  we  consider  the  special  case  o  =  1  + we  find  that  the  ‘'acoustic'  eigenvalue  is 
given  by 

(26)  4  =  ^J{\-qVc^W 

fb'iice,  these  eigenvalues  are  isotropic  in  the  limit  of  M  going  to  zero.  However, 
this  eigenvalue  vanishes  at  the  sonic  line  and  so  the  matrix  is  singular.  In  general,  if 
we  demand  that  the  acoustic  eigenvalues  be  isotropic  then  we  have  a  singularity  at  the 
sonic  line  where  the  eigenvalues  cannot  be  isotropic.  The  two  ways  out  of  this  difficulty 
ar«'  either  to  smooth  the  formulas  near  the  singular  line  or  else  to  giv'e  up  on  isotropy. 

F'or  example  in  [34]  ft  is  chosen  as  zero.  This  results  in  a  ratio  of  about  2.6  between  the 
fastest  and  slowest  wave  speeds  at  A7  =  0.  However,  now  the  formulas  are  regular  at 
the  sonic  line.  This  difficulty  is  not  a  property  of  the  preconditioning  just  presented  but 
applies  ecpially  to  all  preconditioners  e.g.  that  of  Van  Leer  et.  al.  which  will  now  be 
presented. 

The  Van  Leer,  Lee,  Roe  preconditioning  [55]  for  general  non-aligned  flow  in  (^,  du,  di%  dS) 
variables  is 


Pv  = 


-^n/c 


/  M  <  I, 

\  A7  >  1; 

A/  <  1 , 

,  A/  >  1 . 

At  the  sonic  line  /f  =  0  and  r  =  0  and  the  matrix  becomes  singular.  In  both 
thes<'  examples  the  preconditioner  was  constructed  based  on  using  (p,u,,v,S)  as  the 
dependent  variables.  The  rea.son  for  this  choice  is  that  the  matrices  are  essentially 
symiTU’tric  which  this  choice.  However,  if  another  choice  of  variables  is  more  appropriate 
that  introduces  no  difficulties.  Thus,  for  example  [13]  recommend  the  use  of  (p,  u.  i\T) 
variables  for  tin'  Navier-.Stokes  equations.  CJiven  two  sets  of  dependent  variables  w  and 
W  let  VT,,,  be  the  .larobian  matrix  Then,  we  have  dW  =  W,„dn'.  .So  we  can  go 

b('tween  any  sets  of  primitive  variables  or  between  primitive  variables  and  conscTvation 
varialdes.  In  particular  since  the  equations  ar<’  solverl  in  conservation  variables  we  have 
sf'veral  ways  of  going  from  the  primitive  variable  preconditioner  to  a  conservation  variable 


preroiulitioner.  rims,  the  choice  of  variables  used  in  constructing  the  |)reconditioner  is 
dictated  by  mathematical  or  physical  reasoning  and  then  the  preconditioner  can  be 
transformed  to  any  other  set  of  variables. 

•  VVe  can  construct  the  preconditioner  matrix  for  the  conservation  variables.  If  VV 
are  tlu'  conservative  variables  and  w  are  the  primitive  variables  tlu’  Pnni.'ifrradrf  = 

(  H! )  P prnnitirfi  1^  u/)- 

Let  W  <lenote  the  conservative  variables  (p,  m,u, /v)L  with  rn  —  pu.ii  =  pv  .  let  w 
denote  the  primitive  variables  {p,  u,  SY  and  let  w  denote  (p,  i/.e,  7’)h  Then 
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•  VVe  calculate  the  residual  d\V  in  conservative  variables.  VVe  then  transform  dVV  to 
dw  as  before.  Next  we  multiply  by  P  and  finally  transform  back  to  conservative 
variables  flVV  and  update  the  solution.  This  is  algebraically  equivalent  to  the  first 
option  but  requires  three  matrix  multiplies  inst('ad  of  one.  However,  it  offers  more 
flexibility. 

•  Similar  to  the  previous  suggestion  we  calculate  the  residual  dVV'  and  transform 
to  conservative  variables  dw  and  the  midti|)ly  by  P.  At  this  stag<>  we  update  the 
primitive  variables  w.  VVe  then  us*'  the  nonlinear  relations  to  construct  VV  from 
w.  riiis  approach  has  advantages  if  the  boundary  conditions  are  given  in  terms 
of  th(’  primitive  variable's  (p  or  L)  and  so  tln'v  can  be  s|)ecified  exactly  and  not 
approximately. 


These  methods  are  all  equivalent  for  linear  systems  and  tlu'  difference  between  them 
is  mainly  one  of  convenience. 

Based  on  conservative  variables  ('hoi  and  Merkle  [35]  suggest  two  otlu'r  |)recondi- 
tioners.  The  first  is 

1  0  0  0 
0  10  0 

0  0  10 

( M-^  -  1 )  u{M-^  -\)  r( M-^  -  1 )  M-^ 

This  matrix  is  closely  related  to  the  first  preconditioner  Pj  with  o  =  0  after  switching 
between  (p,  u,e,.S')  variables  and  conservative  variables  (see  [51]  for  more  details).  We 
get  a  similar  looking  preconditioner  by  replacing  Ei  in  the  energy  equation  by  and 

tin'll 
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-u  —V  . 

-y- 1  / 

For  the  Navier-.Stokes  equations  they  [13]  suggest  a  different  preconditioner  given  by 

lip  0  0  0  \ 

7^  P  0  ^ 

7i\p  0  p  0 

pu  pe 

( 'hoosing  <^  =  0  or  1  made  very  little  difference  in  their  calculations..  For  inviscid  flows 
■  i  =  .  As  pointed  out  before,  for  both  these  precouditioners  the  ratio  of  eigenvalues  of 
tin'  pri'conditioned  system  is  not  one  in  the  limit  of  M  —  0  but  on  the  other  hand  the 
systems  ari'  not  singular  at  the  sonic  line. 

W'e  thus  again  see  that  the  preconditioner  is  not  unique  for  a  given  set  of  variables. 
Instt'ad  many  matrices  are  capaide  of  reducing  the  spread  of  the  wave  speeds  at  low 
Mach  numbers.  The  main  difference  for  invi.scid  flow  between  all  these  preconditioners 
are  the  eigenvectors  that  r<'sult  from  the  preconditioning.  There  has  been  little  work 
comparing  the  properties  and  efficiencies  of  these  preconditioners. 


3.1  Supersonic  Flow 

VVe  prf'viously  mentioned  that  for  super.sonic  flow  one  can  diagonalize  both  mart  ices 
A  ati  /(  simultaneously  with  a  congruence  transform  (two  dimensions  only).  VVe  now 
('xplicitly  give  this  transformation.  VVe  consider  the  .symmetrizing  variables  (^.  ?/,  v,  P), 
then 
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Let  +  1'^.  We  assume  u  >  0,  e  >  0.  Since  the  flow  is  supersonic  q  >  c.  The 

last  row  and  column  decouple  and  so  vve  consider  only  a  .'3.rd  submatrix.  Define, 
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and  let  Q  =  V^lVi  .  Then 


/  n+ 

0 

0 

(  0  +  -  - 

0 

0\ 

cr-'^Q  = 

0 

u - 

0 

Q*BQ^ 

0  V  - 

uc 

0 

1  0 

0 

tl  ^ 

1  0 

0 

/ 

W('  then  have  the  following  trivial  theorem: 


Theorem  1  If  we  replace  the  matrices  A  and  B  by  the  same  congruent  transformation 
then  this  is  equivalent  to  preconditioning  with  a  non-negative  matrix.  If  the  congruent 
transformation  is  nonsingular  then  the  preconditioning  matrix  is  positive  definite. 

1  lu'  proof  follows  since  QQ*A  =  Q{Q* AQ)Q~'  and  similarly  for  B.  Thus,  the  pre- 
conditioner  P  is  given  by  P  =  QQ*.  The  converse  follows  by  lettiiig  Q  be  the  scpiare 
root  of  P  which  exists  whenever  P  is  positive  definite. 


4  Difference  Equations 

1  ntil  now  the  entire  analysis  has  been  based  on  the  partial  differential  equation.  For  long 
wav('s  it  is  reasonable  to  replace  the  numerical  approximation  by  the  original  differential 
ecpiation.  Since  we  are  mainly  interested  in  wave  speeds  these  are  governed  by  the  low 
frequencies.  It  is  also  possible  to  extend  this  analysis  to  the  finite  difference  approxima- 
ticui.  We  now  make  some  remarks  on  important  points  for  any  numerical  approximation 
of  this  system. 

•  for  an  upwind  difference  scheme  based  on  a  Riemann  solver  this  Riemann  .solver 
shoidd  be  for  the  preconditioned  system  and  not  the  original  scheme.  In  [17]  plots 
are  shown  to  illustrate  the  greatly  improved  accuracy  for  low  Mach  number  flows 
when  the  Riemann  solver  is  based  on  the  preconditioning 

•  For  central  difference  schemes  there  is  a  need  to  add  an  artificial  viscosity.  Accuracy 
is  improved  for  low  Mach  number  flows  if  the  preconditioner  is  applied  only  to  th<' 
physical  convective  and  viscous  terms  but  not  to  the  artificial  viscosity.  Volpe  [bl] 
shows  that  the  accuracy  of  the  origins,  system  deteriorates  as  the  Mach  number  is 


rt'diirecl.  The  aiit  lior  has  had  a  similar  experience  in  thre(' diim'iisional  flews  around 
a  fns(’lage  coidiguration.  Die  use  of  a  matrix  artificial  dissipation  ([a!])  should  hc' 
based  on  llie  pix'conditioiK'd  equations  as  in  the  upwind  difh'renct' scheiiK'.  On  the 
ollu’r  hand  Merkle  (private  commuiucation)  has  indicated  that  he  has  no  difficulties 
with  accuracy  in  tlu'  very  low  Mach  regime.  He  can  take  the  solution  ol)tained  with 
a  prec(jnditionei'  and  us<'  that  as  initial  data  lor  a  nonpr('C(jndit io’-.ed  coch'  which 
tin'll  simply  convi'rges  in  one  time  stej)  with  the  same'  small  n'sidual.  In  this 
cas('  both  the  original  system  and  the  pn'conditioiu'd  system  give  the  same  results 
even  on  the  diffen'iici'  h'vel.  Upwind  schemes  tend  to  have  more  difficulties  with 
accuracy  for  low  Mach  flows  [17]. 

Hence,  Ixith  for  upwind  and  central  differi'iice  sclu'mes  tin'  Hiemann  .solver  or 
artificial  viscosity  should  In'  basi'd  on  p-'IP/t]  and  not  |/1|.  i.e.  in  one  flimension 
solv('  «’(  +  Pfr  =  (|P/'f .  For  a  scalar  artificial  viscosity  |P.'1|  is  rcjilaced 
by  the  spi'ctral  radius  of  P  A  or  <'qui valently  the  time  step  associati’d  with  tin' 
pri'conditioned  matrix.  This  is  erptivalent  to  not  multiplying  the  artificial  viscosity 
by  P. 

•  Similarly,  when  using  characteristics  in  the  boundary  conditions  these  should  be 
bas('d  on  tin'  characteristics  of  the  modified  system  and  not  the  physical  systc'm. 

•  When  using  multigrid  it  is  better  to  transfer  the  residuals  based  on  the  precondi¬ 
tioned  system  to  the  next  grid  since  these  residuals  are  more  balanci'd  than  the 
physical  residuals. 

Preconditioning  is  I'ven  more  important,  when  u.sing  mulligrid  tlian  witli  an  ('xplicit 
schenu'.  With  the  original  .system  the  disparity  of  the eigenvalut's  greatly  affects  tlu' 
smoothing  rates  of  the  slow  components  and  so  slows  down  the  multigrid  mi'tliod. 
[oh], 

•  In  addition  to  convi'rgi'iice  difficulties  then'  are  accuracy  rlifficulies  at  low  Mach 

nundx'i's  [blj.  Sonu'  of  these  can  be  alleviated  by  preconditioning  the  dissipation 
t('rms  as  indicat('d  al)ove.  For  very  small  Mach  numbers  there'  is  also  a  difficulty 
with  roundoff  errors  as  co.  Several  peojrle  have'  suggevsteel  subtracting 

e)ut  a  eeenstant  |)re'ssure'  fnjm  the  dynamic  jrre'ssure.  A  me)re  eletaileel  analysis  [22] 
sugge'sts  re'placing  the  pre'ssure'  ]>  by  p  wlie're'  p  =  J’oVE  f  js  a  repre'sentative 
Mach  nundee'r. 

•  We  concluele'  fre)m  the'  above'  re'inarks  that  the  ste'aely  state  serhitieni  e>f  the  ])ren'on- 
ditiejiie'el  syste'in  may  be'  eliffere'iit  fre)m  that  erf  the  physieal  syste'tn.  Thus,  e)n  the 
fiidte'  elilfe'cence  Ic've’l  the’  pre’ce>tKlitie)iung  ean  improve  the'  ae  curaev  as  well  as  the 
eejiiN'e'rge'nce'  rate'. 


5  Differential  Preconditioners 

In  the'  pre'vions  se'ctie)?is  the'  pre'cejiielit ieuie'r  I’  was  a  matrix.  Feer  the'  neeidine'ar  fluiel 
il\  tiamie-  ('(juat  ions  t  he'  e'le'ine'iit  s  e)f  P  involve'el  t  he*  eie'))i'nele'nt  variable's,  d  here'  are'  se've'i  al 
limitations  with  this  appreeaeli. 


We  first  consider  a  sealar  equation 


(dO) 


u’t  +  (iWj-  +  hwy  =  0, 


We  consider  a  uniform  cartesian  mesh  witli  constant  A.r,  Ai/.  W^e  define  tlx'  aspect  ratio 
for  this  problem  as 

ar  =  aspt'ct  ratio  =  . 

This  can  be  interpreted  as  the  ratio  of  time  for  a  wave  to  traverse  a  mesh  in  the  x 
dirc'ction  relative  to  the  time  in  tlie  y  direction.  W^e  note  that  the  ratio  ^  is  meaningless 
since  this  can  be  changed  by  a  trivial  change  of  variables. 

If  this  aspect  ratio  differs  greatly  from  one  then  the  standard  schemes  will  converge 
slowly  since  a  time  step  appropriate  for  one  direction  is  inappropriate  for  the  other 
direction.  For  a  scalar  equation,  this  is  an  artificial  problem  since,  in  practice,  the  mesh 
would  be  chosen  so  that  the  aspect  ratio  is  close  to  one.  However,  for  a  system  of 
ecpiations  there  are  many  waves.  If  the  aspect  ratio  is  close  to  one  for  one  wave  it  will 
not  be  close  to  one  for  other  waves.  In  the  boundary  layer  for  the  acoustic  wave  ar  = 
('i‘+!-)/Ay  ~  However,  for  the  shear  wave  a**  ~  and  away  from  the  wall  but  in 
tlie  boundary  layer  n  is  much  larger  than  v.  Hence,  any  mesh  that  is  appropriate  for 
th<'  acoustic  wave  is  not  appropriate  for  the  shear  and  entropy  waves  and  vice  versa.  In 
addition  there  are  viscous  effects  that  we  are  ignoring,  so  that  in  practice  the  mesh  is 
constructed  based  on  viscous  effects  and  ignores  both  the  acoustic  and  entropy  waves. 
For  the  scalar  equation  we  are  considering  algebraic  preconditioning  cannot  help  (Li 
and  Van  Leer,  private  communication).  For  a  system  the  preconditionings  we  have 
considered  can  partially  rectify  tlie  difference  of  speeds  between  the  various  waves  but 
do(’s  not  alleviate  the  aspect  ratio  difficulty. 

The  matrix  preconditioners  we  have  considered  until  now  have  a  second  difficulty. 
For  one  dimensional  flow  one  can  choose  the  preconditioner  as  the  absolute  value  of  the 
matrix  A.  Then  all  the  resultant  waves  have  identical  speeds  with  only  differences  in 
the  direction,  positive  or  negative.  However,  in  two  space  dimension  when  the  matrices 
A  and  B  do  not  commute  it  is  not  possible,  in  general  ,  to  equalize  all  the  speeds. 
Fcpiivalently,  we  cannot  diagonalize  the  system  and  reduce  it  to  a  sequence  of  scalar 
equations  even  for  the  frozen  coefficient  problem. 

To  alleviate  these  two  problems  we  shall  allow  the  preconditioner  P  to  contain  deriva¬ 
tive's.  However,  as  before  we  still  demand  that  for  the  symmetric  equations  that  P  be 
invertible  and  be  positive  definite. 

For  the  scalar  equation  (30)  we  consider  a  preconditioner  based  on  residual  smoothing 
[20].  This  is  given  by 


(31  )  (  1  -  6r()rT){  1  -  l^y()yy)RfSy^,r  = 

where  Res  refers  to  the  residual  before  and  after  smoothing.  This  residual  smoothing 
is  usually  introduced  to  improve  the  time  step  and  smoothing  properties  of  an  explicit 
scheme  as  Runge-Kntta  or  Lax-Wendroff.  Here,  we  analyze  the  scheme  from  a  different 
perspective,  that  of  wave  speeds.  We  a.ss>mie  that  the  aspect  ration  for  the  problem  is 
very  large  (i.e.  b  is  large  compared  to  a  or  A?/  is  small  compared  to  A.r  ).  The  question 
we  wish  to  afldress  is  whether  and  i^y  can  be  chosen  so  as  to  reduce  this  aspect  ratio. 
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We  first  consider  residual  smoothing  in  one  space  dimension.  In  this  case  tliere  is 
no  aspect  ratio.  Instead  we  will  show  how  the  concept  of  wave  speeds  explains  the 
phenomena  that  one  should  not  use  residual  smoothing  with  a  very  large  time  step  even 
though  it  can  be  stabilized  by  choosing  an  appropriately  large  3.. 


(1  —  3dj:x)uh  +  nitv  =  0 


We  analyze  this  for  a  semi-discrete  equation  with  time  continuous,  the  first  x  deriva¬ 
tive  approximated  by  a  central  difference  and  the  second  space  derivative  by  a  three 
point  central  difference.  In  order  to  find  the  phase  and  group  velocities  we  consider 
.solutions  of  the  form  w  —  Here  k  is  given  and  we  find  lu  from  the  dispersion 

relation.  For  the  one  dimensional  residual  smoothing  we  have 


k  =  as7n6, 


asiiiO/ \.r 
1  +  23(  1  —  cosO) 


9  =  hAx 

To  find  a  stability  condition  for  a  Runge-Kutta  scheme  in  time  we  maximize  u-’  and 
find  that  the  worst  case  is  co.'?6  =  .  We  then  find  that  the  scheme  is  stable  if 

/^>  j(e-i) 

where  r  =  — .Thus,  from  the  viewpoint  of  stability  we  can  choose  any  time  step 
we  wish  by  choosing  3  sufficiently  large.  Nevertheless,  one  finds  computationally  that 
convergence  to  a  steady  state  is  slowed  down  by  choosing  At,  and  hence  /^,  too  large. 
Optimal  values  are  r  ~  2.  We  shall  now  show  from  the  viewpoint  of  wave  propagation 
that  it  is  not  good  to  choose  a  very  large  time  step. 

Residual  smoothing  adds  a  term  iVni  to  the  original  differential  equation.  Such  a 
term  is  a  dispersive  term  i.e.  the  energy  is  not  reduced  but  now  the  speed  of  a  plane  wave 
is  no  longer  constant  but  instead  depends  on  the  wave  number.  The  main  purpose  of 
this  term  is  to  increase  the  time  stability  limit.  However,  as  in  defining  the  aspect  ratio, 
increasing  the  time  step  is  meaningful  only  if  we  normalize  the  solution  in  some  way, 
otherwise  we  are  merely  rescaling  the  time  dimension.  Hence,  the  appropriate  quantity 
is  not  the  time  step  but  rather  the  time  it  takes  a  wave  to  transverse  one  cell  (assuming 
Ax  is  constant).  The  phase  speed  of  a  plane  w'ave  is  given  by 

u;  a 

^  I  -1-  43smW}2 

Let  3  —  —  1)  and  multiply  Vp  by  r  to  get  the  distance  transversed  in  time  At. 

Then 

2r 

(d2)  .Sp  =  relative  phn.‘<e  di.'itance  =  ^ - - - r - - 

-t-  1 )  —  (r‘‘  —  1  )cc».s0 

For  the  long  wave  lengths  co.sf?  ~  1  and  so  .Sp  ~  r,  i.e.  the  long  wave  lenths  move 
r  times  further  in  one  time  step.  If  we  look  at  0  =  7r/2,  w’e  have  .Sp  =  ^  1-  Thus 

this  frequency  moves  slower  than  witliout  resid)ial  smoothing.  For  the  highest  frequency 
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on  the  mesh  we  have  0  =  ir  and  .s^  =  We  therefore,  see  that  the  higli  freqiienrie<l 
are  actually  slowed  down  by  the  residual  smoothing  and  so  take  longer  to  exit  from 
the  domain,  furthermore  the  larger  At  is  chosen  the  slower  these  waves  go.  Even  more 
important  the  larger  At  is  chosen  the  more  frequencies  that  are  slowed  down  even  though 
the  lowest  frequencies  travel  faster.  The  breakeven  frequency  is  given  by  cosd  = 

We  can  also  consider  the  group  velocity.  For  the  optimal  f:t  this  is  given  by 

do;  +  1  )co.s0  —  —  1  )cos20 

®  dk  [(?-2  -|-  1)  — 

and 

[v^  -I-  1  )cosd  —  —  1  )cos26 

^  ^  ^  ^  {r'^  +  1)  —  —  l)cos0 

The  situation  now  is  even  less  favorable  than  before.  Again,  the  lowest  frequencies  are 
sped  up  by  a  factor  r.  The  frequency  0  =  ir  12  is  slowed  down  by  an  additional  factor  of 
and  the  highest  frequency  0  =  tt  now  reverses  direction  and  goes  upstream  . 

In  figtires  (la-lc)  we  plot  the  phase  and  group  relative  distances  for  r=2,5,10.  As 
demonstrated  above  we  gain  a  factor  of  r  for  the  low  fequencies  but  actually  lose  compared 
with  r=l  for  the  high  frequencies.  As  r  is  increased  more  frequencies  get  slowed  down. 
Because  we  are  considering  the  semi-discrete  equation  and  residual  smoothing  is  purely 
dispersive  there  is  no  damping  of  the  waves.  For  a  Runge-Kutta  scheme  one  finds  that 
as  r  is  increased  that  the  damping  of  high  frequencies  decreases.  Thus,  for  large  r  the 
high  frecpiencies  do  not  propagate  very  fast  and  are  not  damped  either.  This  explains 
one  in  practice  one  chooses  an  r  of  about  two  for  the  greatest  increase  in  the  convergence 
rate  to  a  steady  state. 

We  next  consider  the  two  dimensional  equation.  To  ease  the  derivations  we  shall  con¬ 
sider  the  partial  differential  equation  (dl )  rather  than  the  finite  difference  approximation. 
We  rewrite  (31)  as 


(34)  Ilf  T  f^xf^y^^xxyyt  (IXix  T  ^^^y 

We  are  interested  in  the  effect  of  high  aspect  ratios.  So  we  consider  A.y  <<  A,r  .  By 
rescaling  we  instead  consider  a  uniform  mesh  but  a  «  h.  In  particular  we  shall  choose 
a  =  c,  h  =  \  . 

Consider  solutions  of  the  form  u  =  qj.  equivalently  u  =  f where 

k  =  {kjr,ky)  and  .?  =  (.r.y).  Substituting  this  into  (34)  we  get 


uj{kj.,  ky) 


fkj--\-  ky 

(  1  +  dxki){  1  +  l^yky) 


Hence,  u;(l,0)  =  and  u;(0,  1)  =  If  we  want  these  to  be  equal  then  we  need 

=  0(\),fiy  =  0{j).  This  is  different  than  what  is  normally  chosen  for  in  residual 
smoothing  ([fiO]). 

We  now  consider  differential  preconditioners  for  the  Euler  equations.  We  shall  otdy 
considered  tlie  linearized  ecpiations  with  constant  coefficients.  This  will  now  be  a  matrix 
preconditioner  where  the  elements  of  the  matrix  contain  partial  derivatives.  We  first 

lo 


rewrite  (24)  in  a  more  relevant  differential  form.  Thus,  tlie  Fouler  equations  ran  be 
written  as 

(db)  rt'i  +  hm  =  0 

with  ir  =  {p.u,  V,  Sy  .  We  next  define 

Q  =  udj.  +  vdy 

.Since,  all  coefficients  are  assumed  constant  Q  commutes  with  dj-  and  dy  then 

^  Q  pC^dr  pf^fiy  0  '' 

Q  0  0 

0  g  0 

^  "o  0  a  Q  j 


(.47)  D^Q^-  c^{&l  +  dl). 

We  now  replace  (.45)  by  the  preconditioned  system 

(48)  +  PdLjc  =  0 


/  -pc^drQ  -pc^dyQ  0 

-^-dr  Q^-c^dy  c^drdy  0 


One  can  then  verify  that 


■^dy  c^d:r^y  Q^-c^dj  0 

0  0  0  D 


PDL  =  gDI,  Po'  =  D-'Q-'L 


One  can  of  course  replace  the  D  in  the  lower  right  corner  of  Po  by  the  identity 
matrix.  Then  PdL  is  not  the  identity  matrix  but  is  still  a  diagonal  matrix.  We  can  use 
simpler  matrices  than  Pp  by  considering  congruent  transformations.  We  consider  the 
symmetrizing  variables  c(^,  «,  r,  S),  theti 


Q  rdj.  cdy  0 
cd,  goo 
cdp  0  g  0 
0  0  0  g 


-r()y  0 

0  0 

1  0 

0  1 


g  0  0  0 

-ctK  1  0  0 

—  cdy  0  1  0 

0  0  0  1 
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/  DQ  0  0  0  \ 

p^r  pt  ^  0  0  0 

0  0  0 

\  0  0  0  g  / 

so  w('  liavo  diagoiializf'd  L  l)y  a  rongriu'nt  transformation.  Hut. 

PeLP^  =  P|:‘'(P^PeL)Pe. 

so  tli('  roiigriu'iit  transfortn  is  similar  to  a  precoiulitioning  with  a  positive  dt'finite'  matrix 
PePe  Alt(>rnativ('Iy,  (  P^PeIL  is  similar  to  a  (liagonal  matrix. 


PePe 


(  Q  -c()y  0 

_  —c()j.  1  +  r^()j  c^dj-dy  0 

-('()y  C^()j()y  1  +  C^dl  0 


Notr  that  PePe  looks  similar  to  Pd  Imt  is  not  idiuitiral.  PePe  has  fewer  deriva¬ 
tives  along  the  ideiitieal  hut  P^PeL  is  only  similar  to  a  diagonal  matrix  while  PqL  is 
diagi)nal  and  evaui  a  scalar  differential  operator  imdti])lying  the  identity  matrix.  These 
transformations  are  indep('ndent  of  the  flow  regime  as  long  as  the  preronditioner  is  non¬ 
singular. 

T1h'S('  |)reconditioners  are  ('onn('ct(>d  with  the  techniqiu's  used  in  distributive  (lauss- 
.Seidel  smoothers  for  nndtigrid  methods  ([f)|4"])- 

It  remains  to  show  that  P  is  nonsingidar.  VVe  have  four  eigenvalues  and  corresponding 
eig('nfunctions.  As  usual  tluMuitropy  wave  decouph's.  For  this  wave  P  has  an  eigenvalue 
I)  and  an  eigenfunction  (0,0,0. 1)  .  For  the  shear  wave  P  has  an  ('igt'nvalue  1)  and  the 
('igenvector  is  (ei,r2.e3,0)  where 


Di'\  =  0 

dv) 

(fX  Ojl 

The  other  two  ‘acoustic’  eigenvalues  of  P  are  ±  rC 
satisfy  the  pseudo-difh'rential  erpiation 


)f.  Oy  and  the  eigenvectors 


We  therefore  have  to  show  that  the  eigenvalues  ar<'  all  nonzero  so  that  P  is  nonsin¬ 
gular.  file  opc'rator  1)  is  just  the  i)otential  operator  i.e.  for  any  variahh'  w 


Dll'  =  (ir  —  c^fu’rr  +  2uru\„  +  (v^  —  r^)ir 


For  siihsoiiir  flow  litis  is  an  ('Ilipt ic  i)p<'ralor  and  so  invortible.  For  siiporsoiiic  flow  I) 
is  a  hyix'rholic  ojx'valov  .  Similarly,  Q  is  a  hyporbolic  operator  di'iiotiiig  ronvoction  along 
a  streainliiK'.  d  ims  gi\('n  appropriate  boundary  eoiiflitions  this  loo  slii/nlfl  be  inva'rt il'<l(\ 
At  a  stagnation  point  Q  is  singular  and  so  it  is  nee<'ssary  to  limit  the  vabn's  of  n  and  in 
V  in  the  definition  of  Q  so  that  tlnw  do  not  beeome  loo  small  in  a  neighborhood  of  the 
stagnation  point.  A  similar  smoothing  is  needt'd  near  the  sonic  line.  These  arguments 
have  Ix'tMi  applied  to  Pd  but  sindlar  arguments  work  for  P^Pe- 

W  ith  residual  smoothing  and  Pq  or  PePe  ''■<'  have  increased  the  oiah'r  (T  the  system 
and  so  chang<'d  the  nnmlx'r  of  boundary  coiuiitions  needed  for  the  equation  to  Ix'  well 
pos('d.  'ho  avoid  this  difficnlty  w<'  do  not  solve  tlu'  equation  (38).  Instead  these  })recon- 
ditioners  are  nst'd  as  a  post  processor  for  tin'  nsnal  Knler  or  Navier-Stokes  erpiations. 
rims,  at  ('ach  time  st('p  wv  calcnlate  a  residual  based  on  one's  favorite  schem*'.  This 
giv('s  a  pr('(licl('d  value  of  the  change  in  liiiu'.  Afe,,r,,/„(,vy-  Wt'  also  update  I  lu'  boundary 
conditions  for  tin'  standard  fluid  dynamic  ecinatiinis.  W('  then  operate  on  Aic  with  P 
with  the  boundary  condition  that  AnVorrrOff  =  0  -  >  e.  we  don't  change  the  boundary 
valiK’s  calculated  by  the  predict t)r.  W  hen  we  reach  a  st<’ady  state  for  the  fluids  equations 
we  an'  solving  PAn’,.„rr,,f,./  =  0  with  Z('ro  boundary  conditions.  Since  P  is  invertible 
Au’  =  0,  i.e.  we  pres('rve  the  steady  state.  Thus,  in  ('sscmce  we  are  imposing  the  fluid 
dynamic  l)oundary  conditions  betwcxm  the  P  operator  and  the  L  ojrerator. 

6  Alternate  Methods,  Time  Dependent  Problems 
and  Viscous  Problems 

rile  justification  for  preconditioned  schemes  began  with  low  Mach  number  flows.  For 
such  flows  other  t<'chni(iues  exist  beside  preconditioning  tlie  ecjnations.  The  method  of 
tilin'  inclining  has  similaritii's  to  preconditioning  [loj  . 

rhe  basis  of  one  such  method  is  to  u.se  an  iiiqilicit  scheme.  Ho.v’ever,  a  two  dimen¬ 
sional  implicit  method  is  too  «'X[)ensiv<'  to  lie  efficient.  Thus,  one  classically  uses  an 
.\f)I  approach.  llowc'Vi'r,  it  is  known  tliat  with  .ADI  one  cannot  choose  a  very  large 
time  st('|)  and  convi'rge  (piickly  to  tin*  steaily  state.  The  sjilitting  errors  that  occur  in 
the  ADI  method  couples  the  wavi's  together  aiul  one  cannot  choose  an  appropriate  time 
step  for  ('ach  wav('.  Instead  one  attempts  to  separate  thos('  terms  in  the  equations  that 
contribute  to  the  fast  acoustic  wave's  from  the  slow  components.  One  than  can  use  a 
s('mi-im|)licit  method  which  is  inqilicit  for  the  fast  waves  and  explicit  for  the  slow  waves, 
riins,  the'  stability  limit  of  the  scheme  is  goveriu'd  by  the  convective  speed  rather  than 
the  acoustic  speed  .  Tin'  ('xplicit  part  can  be  ('ither  a  leapfrog  method  ([20],  [21]),  or 
a  two  st('p  iTK'thod  [22],  This  can  also  be  exte'iidi'd  to  the  Navie'r-Stokes  equations  [23], 
.\lt('rnal iv('lv.  onc('  tlu'se  components  are  identified,  one  can  split  the  ('qnations  in  several 
|)i('ces  and  solv('  ('ach  one  s('|iarately  as  in  tin'  classical  splitting  nu'thods  [2]  .  In  this 
cas('  OIK'  can  ns('  an  implicit  iin'thod  for  tin'  fast  waves  and  an  explicit  method  for  the 
slow  waves  and  in  addition  oik'  can  sjilit  olf  the  viscous  t('rms.  d'lu'.se  methods  work  for 
both  t ime  depi'inh'iit  and  sti'ady  stat('  probh'iiis. 

.A  difb'n'iil  alternative'  is  to  add  l('rms  to  the  ('qualions  which  disappear  in  the  sti'ady 
state.  This  has  a  coniK'dion  with  pri'condif ioni'd  iiK'thods  wln'ii  time  derivatives  are 
add('d  to  the  ef|uations.  Ilowevi'r.  in  this  approach  other  terms  can  be  added  bt'side 
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tinu'  (Ifiivat i\('s.  One  ('xampu',  is  ft)  assiiiiK'  that  tlic  total  ('iithal|)y  is  constant  in  tin' 
steady  stat(‘  for  t  lu'  comi)r('ssil)l(‘  iiniscid  r(|nations.  Oiu'  ran  tluMi  add  terms  to  the 
('f)nations  that  depend  on  the  df'viation  of  tlu'  rurr<Mit  <'ntha!py  at  ('arh  point  from  tliis 
constant  st('ady  statt'  entlialpy  [2")].  h'or  tlu'  incompn'ssil)l(' ('(piat  ions  one  can  add  tlie 
di V('i t>;ence  of  t h(‘  \'elocity  fitdd  or  time  diM'ivat i ves  t)f  tlie  dixangt'iice  to  the  motnentnm 
eipiation  [11].  [Id]  .  Oih'  can  also  consid<M'  a  more  ,u;enertd  eipialion  ol  state  that  rediicc's 
to  the  pliNsical  oik'  at  the  steady  state*  [l  1]  .  In  [’27]  they  atialyze*  the*  General  case  of  such 
pseinlo-mist eady  s\st ems. 

An  extension  of  tins  techinepie  is  to  modify  tin*  dilh'rent iai  eepiation  to  remove*  the* 
aeenistie-  \vave*s  eer  e)the*r  ’hael'  fe*atnre*s.  One*  must  th<*n  justify  that  the*  seelutieens  eeh- 
taine*el  te)  the*se*  me)elifie*d  e*e|uatie)ns  are*  edeese*  tee  the*  oriftiind  e*epiatie)ns  for  seeme*  flow 
re*(>,ime'.  l  ypie-al  e*xan\ple*s  are*  t  he*  varieeiis  I.eew  Mae  h  numlie*r  e*xi)ansions  feet  t  he*  Iluiel  ely- 
namie-  e^epiatieuis  or  the*  f;;e*e)st  rojehie-  e*epiatie)ns  as  an  aiepreeximat  ieen  tee  the*  shalleew  wate*r 
e-epiatieens  in  tne*t e*e)re)le)Ky. 

h'eer  ine'otnpre*ssil)le*  lleew  peepidar  se'he'me*s  are*  the*  SIMPLI''  [.dfl]  tuiel  MAC  [IS]  algee- 
rit  hms  and  t  he*ir  p;e*ne*rtdizat  ieens.  rhe*se'  usually  re*fpiire*  t  he*  seelut  ieeii  eef  a  Peeisseen  e'epiat  ieeii 
fejf  the*  pressure*  anel  the*n  a  pre*ssur(*  ee)rre*e't ieeii  is  use*el  tee  u])elate*  the*  me)me*ntum  e'epia- 
t  ieens.  rhe*se*  me*t  heeels  e-an  t  he*n  he*  <j;e*ne*ralize*el  te)  t  he*  e  e)mi)re*ssil)le*  e*epiat  ieens  [21].  Me-rkle*. 
Ve'nkate*swaran  ain!  Hue*le)W  [.'{7]  ee)mpare  sue  h  me-theiels  tee  the*  pre*ee)nelit ie)ne*d  te*e  hid<pie*s 
elise’usse*el  in  this  |)ape*r.  We*  agaiti  st  re*ss  that  the*  elilfe're'nee*  in  t  he*se*  appre)ae’he*s  is  ne)t 
whe'the*r  ele*nsity  or  i)re*ssure'  are*  use'el  as  the*  el(’pe*nele'nt  \;irial)le'  as  eene*  ean  trtmsfortn 
he*twe'e*n  the  se*  variahle-s.  d  ims,  feer  e'xtimple*.  eeiie*  ean  moelifv  the*  eeempre'ssihle*  eont innity 
e*epiatie)n  hy  rcplaeinc;  the*  time*  ete*rivat  ive*  of  the*  e|e*nsity  with  a  time*  de’iivat  ive*  e)f  the* 
(U’e-ssure*.  I  his  is  just  ane)t)ie*r  e*\ample*  eef  a  matrix  pre*ce)nelit ie)niti,e,  as  eene*  e  an  e*xpr(’ss 
the*  pre*ssnre*  e|e*rivat i ve*  <ts  a  eeemhinat ieeii  eif  a  e|(*nsity  eh'rivat i ve*  te)j>;the'r  with  me)me'ntnm 
anel  e*ne*r,2;y  de*rivat ive-s.  .As  e|e'se  ril)e*e|  aheive*.  it  is  a  prosi'tiinmins  de'e  ision  whe*the*r  one* 
sheeulel  use*  this  inodificel  e*eiuation  te)  u|)elate*  the*  i)re*ssnre' ainl  t  he*n  transfe*r  te)  eh'tisity  or 
to  e  ale  nlate*  the*  the*  ap|)re)priat e*  pre*ee)ndit ieening  matrix  ttttd  u|)elate*  the*  eh'nsity.  I*'e)r  a 
litie-ar  syste-m  the*  two  appre)aehe*s  are*  iele*nt ieal. 

For  time*  e|e*pe'nele'tit  prol>le*ms  the*  first  appre)ach  just  dise  usse*el  is  use*lid.  lle)We've'r. 
the*  pre*ee)nelit  ieeiie’el  me’thoels  and  the*  se*ee)nel  appre)aeh  e)f  this  se*etie)n  ele*stre)y  the*  titne* 
ace  iiraey  utih'ss  the*  coe'ffie  ie*nt s  of  the*  pe-rt urhat ieen  are*  che)S('n  as  a  functie)n  e)f  the’  tne*sh 
size*  anel  se)  otdy  affe-ct  te-rtns  of  the*  e)rele*r  e)f  the*  ae  e  nrae  v  of  the*  se  he*me*.  .A  tnore*  popular 
appreeaeh  has  l)e*e'n  te)  use*  a  twee-lime*  sehe*me*.  In  this  appreeaeh  e*arh  ne*w  time*  le*ve*l 
is  ee)nsiele*re*el  as  the*  ste*aely  stale*  of  seeme*  |)re)l)le*m.  .Alle*rnal  ive*ly,  the*  i)hysieal  time* 
ele*ri\at i ve*s  are*  ee)nsiele're*el  a  forcin.e;  te*rms.  One*  neew  nse*s  the*  pre'eeenelit ie)ne*e|  Tne*thoels  tee 
ae  hie'Ve*  this  'ste'aely  state*'  whie  h  in  re*ality  is  the*  se)hilie)n  at  the*  ne*xl  time*  ste'p.  Ile*nee*. 
the*re*  is  the*  lehysieal  time*  I  anel  an  artifieial  time*  r  anel  r  p;e)e*s  lee  infinity  as  an  inne-r 
le)op  within  e*aeh  time*  ste*p.  [12]  ,  [17],  [h'']).  I  hus, 


ih  ^  m 


()  f  ()(/ 

—  -f  —  =  0. 

().r  ()ti 


I  he*  main  eliflieidty  with  this  appreeaeh  is  its  (*ffieie*ney.  It  is  re*ase)nal)le’  te)  use*  sueh 
a  te’elnnepie*  oidy  if  e*ach  'ste*aely  state*'  pre)l)le*m  ean  he*  se)lve*e|  with  little*  e*ire)rt.  One* 
ael vant a,e;e*  is  t hat  one  usually  has  geeeeel  initial  ,e;ue*ss  feer  t he*  seelnl  ieeii  hase-d  een  t  he*  seel ut  ion 
at  pre’\ie)us  time*  ste-ies.  Ile)we*ve*r.  it  typieally  take*s  10  suhite'rat ieens  feer  e-aeh  lime*  ste*)). 


Heine,  this  approai  h  is  ten  times  inort' expensivt'  lliaii  a  straiglit  implicit  method.  One 
can  also  use  an  .Xmvton  itmation  [dS]  at  each  tiiiu'step.  nevertheless,  a  stmii-itni^lieit 
appi'oach  in  ([20]  -  [2d])  seems  at  t  raet  i\-('. 

All  the  methods  diseiissc'd  thus  far  hav<'  het-n  has<'d  on  an  inviseid  atialysis.  For 
the  .\a\ ier-Stoki's  e(piations  at  hi,i>;h  l^eynolds  numlx'r  \v(’  do  not  expc'ct  atiy  itnportant 
ehanyes  outside  t  lu'  houinlary  layi-r.  Insiih'  the  l)oimdar\’  layc'r  r  isetnis  eih'ct.s  modify 
the  ei[ieu\alues  of  the  dilh'fent  ial  opi'rator.  \V<'  t  hu.s  wish  to  (xpialize  t  lu'  eoiit  rihiit  ioti 
of  three  (|uantities.  the  aeoustie  wa\'es.  the  eoii\eetive  waves  and  the  riseous  tertns.  In 
particular  the  \iseous  eigetivalues  an'  vc'iv  stiff  and  so  the  eigenvaliK's  of  the  solution 
operator  are  no  longer  well  eoinlit  ioiu'd.  .All  the  preeoinlil  ioiu'rs  pr('S('tited  above  d('pend 
on  free  paramett'is  (d.o.r.^)  .  ()|)timal  value's  for  the'se'  parattK'tc'rs  we're  give'ii  feer 
invise  id  flow.  ,\  sim|)le’ exte'iision  of  the'  alteeve'  me'theeels  tee  vise-e)us  fle>w  weiulel  ke'e'p  the' 
same'  form  for  the  pre'e  ondit  ioning  matriee's  hut  alleew  the'se*  paratne'te'is  tei  alsei  elepe'iiel 
on  the-  Ke'v  Holds  of  Fraud  1 1  numhe'r  (se-e-  for  e'xami)le'  [  10]  .[  hd]  ).  dduis.  feir  e'xample'  e)ne' 
finds  that  for  the'  eiriginal  pse'inlo cotnpre'ssileility  me'theeel  that  .i  she)nlel  ine're'a.se'  as  the 
i{e'ynolds  numhe't  is  eh'cre'ase'd.  In  [Id]  a  ne'w  pre'e  ejtielit  ie)ne'r  is  int  renlue'e'el  Mase'il  on  a 
lihysieal  analysis  eif  the'  .Navie'i-Steike's  e'epiations  (se'e'  2!)).  1  he'  eliflie  ulty  is  that  the'  time 
ste'|)s  are'  gove'riH'el  hy  the'  aeeiustie  anel  eon  ve'et  i  ve' spe-e'els  atiel  alse)  a  vise'e)us  e'e)ittril)utie)n. 

.\  basic  preible'tti  feer  the'  pre’conelit  ieene'el  .X'avie'i-.St  eeke's  e-epiat  ions  is  we'll-peeseelne'ss.  Fe)r 
t  he'  in\  ise  iel  e'epiat  ieins  one'  e  an  she>w  t  hat  wit  li  t  lie'  pre'e'onelit  ieine'r  Pj  '  bat  the  e'fpiat  ions 
(an  lie-  symtiK't rize'el  if  <\..i  satisfy  the'  ine'ejuality  (lb),  (se'e'  [.">}]  ).  d'lie'  pre'reinelitieine'r 
Pv  is  const ruete’d  from  the  symittel fie-  feirm.  Ile-nee',  in  beith  ease's  we-  rati  re'eluee  the' 
pre'condit ione'd  eepiations  to  a  symtne'trie-  hype'ibeilir  .syste'tn  atid  se)  it  is  we'll  pose’el.  Otice 
oti('  aelds  the  viscous  te'ritis  this  atialysis  is  im  leinge'f  valiel.  One*  peissibilit y  is  to  start 
with  a  form  that  is  symtne'trie  feir  beith  the'  iiiviseiei  atiei  viscenis  te'rms  [1],  If  otie'  uses 
a  |)ositive  de'finite'  pree ondit ione'r  for  the'se'  variable's  the'ti  statidarel  e'lie't'gy  argutiu'tits 
shows  that  the'  line'arized  jire'cemdit iotte'd  syste'iii  is  we'll-peise’el. 

We  neiw  analyze-  the-  pre'e  emdit ione-r  P;  a  little  meire-  carefully  for  the-  inconipressilde 
.Navie-r-Stokes  e(|uations.  We-  also  line-arize-e'  anel  so  the-  e'eie'lfie  ie'nts  i/.  i\  d  are  eotisiele-re-d 
as  constant.  I  he'  re'sultant  jire-eondit ieine-el  e-epiat ie>ns  are 

~JjPi  +  "<■  +  '’v  =  b 

oil.,  . 

(11)  ^  ~ 

oi'o  ,  ,  . 

"b  l‘ll''r  +  '’(iCy  +  />„  =  //Wc 

\\('  tie’xt  diffe'ie'tit  iate' t  he’  se'coiiel  e-epiat  ie>n  wit  h  re's|)e'ct  t o  x  and  t  lie- 1  hirel  wit  h  re'spect 
to  y.  We- re-place- t  he’ dive'ige'iie  e' of  t  he- ve’leicity  fremi  t  he- first  e-ejuatioti.  Let  H  =  t/o(  t/j-r + 
c,„)  +  C|i(";v  +  I'’!!/)-  I  he-n  the-  lue-ssure-  p  satisfie’s  an  aeenist iedike-  e-epiatiein 

(  12)  - 77/'"  “b  d"  7T7("<)/».,  +  f’o/'v)/  =  —H 

riius.  we-  re-|)lace-  the-  Feiissein  e-epiatiein  use-el  in  the-  MAC  type-  approach  fry  a  ge-ne-r- 
alize-d  wave-  eepiation  for  the-  pre-ssure-.  We-  Feuirie’r  transfeirtn  (12)  ,  i.e-.  p  —  , 
and  |F|^  —  ky  +  k\.  I  he-n. 


(i;{) 


VVf  first  coiisich'r  tlit'  ease'  o 
primitive'  ('filiations).  TIk'h 


—  |fv(i/nA'i  +  I’oA'i)  +  ^  —  ff- 

0  (i.f'.  t lu'  original  psoudo-compn'ssihilit y  for  the' 


(  M)  ou'  —  - - - 

We  now  liav('  two  rf'ginu's  to  considf'r 

case'  1:  \k\  small  (i.('.  \k-\^  < 

d  lu'ii  (41)  givf's  u,’  .  As  ('xpe'ctf'd  //  int rodiKX's  a  di'cay  in  tlu'  aroiistir  wave'  .  1  lu' 

s|)('('d  of  tli('  wave'  (re^al  part  of  u,’)  is  now  slowe'd  down  for  the'  same'  ,i.  Wv  thus  shoiiM 
choose'  a  large'r  i  as  /i  incre'ase's  to  compe'iisate'  for  this  (see'  also  [l-l]). 

case'  [A'l  large'  (i.e'.  |A'|^  >  Now,  u;  =  ///|A'i^  Jl  ±  ^  1  —  ld'^//‘^|A’|'‘^j .  lle'nce', 

u.'  is  pure'  imaginary,  d  ims,  tlu'se'  high  fre'eiue'iicie's  do  not  propagate'  and  their  damping 
is  re'duce'd  hy  ,4  (for  the-  smaller  damping  mode').  1  hus.  one'  also  wants  to  incre'ase  .i  so 
that  most  of  the'  mode's  in  the'  domain  corre'spond  to  small  |A|  . 

\V('  iK'xl  eon.side'f  non-ze'ro  o.  Le't  y  —  +  coA’i) 

y  +  //<|A'!^  ±  -  //'^lA'I'  +  ‘i/y/zlAi-' 

(do)  i - - - . 

d’aking  real  and  imaginary  [larts  of  the  scfiuirr  root  we'  se'e'  that  only  y^  I'lite'is  into  the 
imaginary  part  of  u.',  i.e'.  the'  de'cay  rate'.  So  the'  sign  of  n  is  not  important  for  viscous 
('[fe'cts.  dhiis,  it  se'ems  that  o  has  no  major  impact  on  viscous  flows  and  its  advantage 
come's  from  e(|ualizing  the  flow  spe'e'ds  of  the'  inviscid  iiortion  of  the'  flow. 

7  Computational  Results  and  Conclusions 

.XuMU'rous  authors  have  use'd  some' of  llu'se'  preconditioners  for  lioth  incomiue'ssihh'  and 
compressihle  flows.  A  se'le'ction  of  |>ap('rs  is  pre'se'uted  in  the'  hililiography.  Ile'ie'  we' 
summarize  a  few  of  tlu’se  calculations.  Most  of  the'se'  computations  have'  use'd  ce'utral 
diffe're'iice'  approximations  of  the  spatial  de'rivative's  and  e'ithe'r  a  Huuge-Kutta  e'xplicit 
scheme  or  an  A.*^).!.  implicit  sche'iue  in  time. 

For  the  original  pseudo-compre'ssihility  ('e|uations  a  numlx'r  of  authors  (e'.g.  [10], 

[to],  [11]  )  have  found  that  a  constant  ,4  works  Ix'st,  Hizzi  and  Kriks'^on  [fy]  suggest 
=  ?uu.r(0.d,  r( +  e^))  with  1  <  r  <  5,  se'e'  also  [9]  .  lu  [-48]  the'V  also  ('X|)lor('  similar 
issue's  with  re'gard  to  upwind  scheme's.  As  Ix'fore'  tlx'ir  constant  O.tl  must  de'pe'ud  on 
the'  normalizations  use’d.  Arnone  ([4],  [4])  has  use'd  the  original  pseudo-compre'ssihility 
me'thod  to  solve  inviscid  and  viscous  iiie-ompre'ssihle  flow  about  cascade's.  .A  Hunge'- 
Kutta  method  is  used  which  is  accele’raled  by  a  multigrid  te'e  hniepie'.  Fins  me'thod  has 
be'e'ti  e'xte'ude'd  by  the'  author  to  inelude'  the'  pre'conditione'r  Pj.  In  the'se'  calculations 
W('  find  that  ,i  =  constant  is  more’  robust  than  (hex)sing  .i  to  de'pe'iid  on  the  spe’e’d 
of  the’  flenv.  In  meest  case’s  using  a  variable'  .4  eause's  the’  iterations  to  dive’rge'  though 
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wIk'ii  tlu’V  iK)  <()iiV(M(>;f'  it  is  faster  than  the  constant  .i.  I’aul  and  ('arlsi)n  [10]  have  a 
similar  three  dinu'iisional  code  for  exti'inal  How  over  win,u;s.  I  his  code  has  also  hec'ii 
exteinh'd  to  inchnh'  Pt  In  both  these  codi's  the  conv('rG,enc('  is  also  vt'iv  dependent  on 
tin'  houndary  conditions  imposed,  lor  sorin'  honndary  coinhtions  the  code  con\’eri2;ed  for 
a  ranye  of  a  ami  tln'ii  a  =  I  gave  the  fastest  converg('nc<>  rate's  as  ('xpe'cti'd.  Ilowe've'r. 
for  othe:-  homnlarv  ceinditions  only  tlu'  original  pse'inki-compri'ssihility  method  a  =  0 
would  conV('rg('.  It  is  suspecte'd  that  t  lu'  dillicult  ie's  are'  eonne'et  e'll  with  init  iaiizat  ieui. 

1  hus.  a  =  I  though  faste'r  ma\'  In'  h'ss  rerhust  .  It  weiiilel  t  he're'ieiie'  he'  ne'iee'ssary  tei  start 
the-  calculation  with  a  =  0  anel  only  onee'  the-  asympteitie-  re'gieui  is  re-ache'd  to  change'  to 
a  =  I . 

Ilsu  [10]  also  salve's  t  In' ini'ompre'ssihle' e'epiat  ions  using  Px-  In  t  his  case' an  u|)winde'd 
ai)|)roximat ion  is  use'd  and  tin'  seilution  is  aehanee'el  using  an  .A.D.l.  me'theiel.  Tln'v 
e'xamine'  in  more  de'tail  the'  mihienee'  e»l  e>  anel  i.  Due'  tei  the'ii'  implie  it  seilve'r  the'  eoeh' 
eeem  e'fge'uee'S  ill  all  the'  ease's  t  he'V  trie'll,  mainly  Hows  about  a  ili'lta  wing..  Howi've'r.  t  he'y 
also  liinl  that  1  is  faste'r  than  the-  variable'  -i.  'I'lie'V  prineipally  inve'st igate'el  ei  =  —  1 

but  iinlii  ati'  that  othi'r  a's  be'ha\e'el  similarly.  I’ln'ie'  have'  be'e'ii  no  eeunputat ions,  tei  elate', 
for  t  In'  ini'oinpre'ssible'  e'epiat  ieiiis  using  t  he'  Py  pre'e'onelit  ieine'r  iliie'  to  t  In'  ne'wne'ss  eif  t  his 
apiu'oai'h. 

bar  the'  |■ompre'ssible' e'epiat  ions  at  leiw  .Mai'h  numbe'is  e'arly  ealciilat  ions  we're'  eleuie' 
by  Hrih'v.  .Mi  Deinahl  ainl  .Shamroth  [S]  anel  a  late'r  by  D.  Choi  ami  .Me'rkh'  [11],  anel 
also  ^■.ll.  Choi  anel  .Me'i  kle'  [H  i]  .  d  he'se'  me'llienls  have'  mainly  use'll  .A.D.l.  nu't  hoels 
though  some'  re'sults  with  Hungi'-Kutta  se  he'ine’s  have' alsei  be'e'ii  ae'hii've'il.  .Meire'  rei'cutly 
([17], [■)")])  ri'siilts  have'  be'e'ii  aehie've'il  with  the'  Py  pre'e'oinlil  ioiii'f  in  e-onjune-t  ion  with  an 
u))\vinil  si  lu'ine'.  (ioilfri'V  (private'  e  eunmunie  alion)  inilie  ate's  tliat  tlii're'  is  not  a  gre'at 
iliire'ie'iice'  be'twe'e'ii  the'  twei  pre'e'onelit ioiie'rs.  The  u.se'  e)f  the'  eeirre'ct  Hii'tnann  solve'r  was 
more-  impeirtant  than  the’  ele-tails  of  the'  pre'e-einelit ione'r. 

.\Iuih  of  the'  most  rei'i’iit  work  has  geme'  inte)  e'xle'iieling  t  he'si'  ri’siilts  to  the'  Na\ie'r- 
Stoke's  e'epiations  [IH]  ami  ehemistry  ([IT],  [iSj.  [AS]).  .A  number  of  authors  have'  also 
invi'st  igati'il  e’xte'iisions  tei  t  inn' ele'pe'iiele'iil  proble'ins  base'el  on  a  t  wei-t  inn' aiipreiaeh  ([lb], 
[IS],  [(;-]). 

lie'!!'  wi'  pre'se'iit  eiiily  one'  se't  of  re'sults.  lids  is  feir  ini'om|)re'ssible'  How  aremml  a 
\  Kl  1  asiaele'  with  a  non pe'iioilii'  ine’sh  aeross  the-  wake'.  1  he'  me'sh  is  shenvn  in  figure' 
■Ja.  .A  Huiige’  Kut  ta  multistage'  seln'ine'  is  use-el  with  a  multigriel  ae'ee'le'iat  ieui.  1  In'  I'eieh’ 
is  a  e'Xte'iisioii  of  that  eif  .Arne)in'  ami  Ste'e  ee)  [ij.  I  In-  He)W  is  tiirbuli'iil  with  a  He'N  iiolils 
numbe'r  eif  ."lOO, ()()()  anel  Halel win- Lomax  type-  lurbule'iice'  moili'l  is  use'll.  In  table'  I  we' 
pri'si’iit  the'  ri'siiiual  of  the’  jire'ssiire'  afte-r  -AO  sle'ps  on  the'  lirst  me'sli.  •)()  ste'|)s  eiii  the' 
si’e’oml  mi'sh  ami  AOO  ste’ps  on  the'  line'st  nie'sli.  We'  thus  se'e'  that  n  =  1  gave'  the'  faste'st 
e'on\ I'lge'in  1'  rate's,  though  the'  eliire’ri'iiee's  we’ie'  not  ve'iv  large'.  We'  we're  able'  tei  run  only 
the-  moilitie'il  \'an  Le-i'i  e't.  al.  lue'ieimlit ioni'r  ami  e've-n  that  only  with  a  emnstant  f>  anel 
.i  wit  h  e>  =  [  as  eippose’il  tei  t  he'  \alue'  e)l  n  give'll  in  ('d'd).  W  it  h  t  his  \alue'  ol  o  t  In'  te'i  ins 
with  1/“’  -f  c"’  ill)  not  ap})i'ar. 


method 

S 

o 

residual 

precoiulition  AV 

1 

1 

(i.tiij 

precondition  AV 

1 

0 

6.07 

precondition  .AV 

1 

-1 

5.76 

no  precondition  AV’ 

1 

1 

6.  .44 

no  precondition  AV 

1 

0 

6.24 

no  precondition  A\^ 

1 

-1 

5.76 

precon  AV  t'tp  (22) 

1 

0..") 

6.44 

no  precon  .AV  erj.  (22) 

1 

0.7) 

6.22 

Tal)lo  1:  Convorgonrt' rato 


III  figure  (2e)  we  also  plot  tlie  roiivergeuce  rate  for  the  first  example  in  the  table. 

Ill  eoiielusion  these  eompiitatioiis  show  that  one  ran  ralnilati'  both  iiivisrid  and 
visrons  flows  and  evani  those  with  rliemiral  reactions  over  a  large  range  of  Mach  numbers 
going  down  to  A7  =  10'"’  in  some  cases.  There  is  need  for  further  work  on  the  effect 
of  till'  parameters  in  the  preconditioners  on  ihe  convergence  rates.  It  is  not  understood 
why  constant  li  seems  to  ire  the  fjest  choice.  There  is  also  need  for  further  investigation 
on  the  effect  of  boundary  conditions  on  the.se  preconditioners. 
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Figure  2a:  Mesh  for  VKI  cascade 
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Figure  2c:  a=1  p=1  Re=500,000  precondition  viscos 
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